Relationship between tectonic tremors and 3-D distributions of thermal structure and dehydration in the Alaska subduction zone

The Alaska subduction zone is characterized by a subducting oceanic plateau, which is referred to as the Yakutat terrane. Tectonic tremors occur in this zone, and there are few volcanoes above the subducted Yakutat terrane. In this study, we performed a 3-D numerical simulation of a thermal structure associated with the simultaneous subduction of the Yakutat terrane and Pacific plate to elucidate the mechanism of tectonic tremors, which typically involve the presence of water. We calculated the water content distribution near the slab surface by using the thermal structure obtained from our simulation and phase diagrams of the hydrous minerals included in the slab. As a result, dehydration from the marine sedimentary layer and oceanic crust was observed near the area where tectonic tremors occurred. Tectonic tremors occur only in the Yakutat terrane because the marine sedimentary layer and oceanic crust are thicker there, and the total amount of water content in these layers is higher; therefore, the amount of dehydration is also higher there than in the Pacific plate. Additionally, there are few volcanoes above the subducted Yakutat terrane because little water remains within the slab beneath the volcanic chain where magma is produced.


Text S1
Model Fig. S5 shows the 3-D parallelepiped model domain. The +x axis is taken to be perpendicular to the trench axis in the direction from inland to the trench, the +y axis is oriented along the trench axis from the Yakutat terrane side to the Pacific plate side, and the +z axis is vertical in a downwards direction. The model size is 600 km × 700 km × 200 km (X × Y × Z) and the number of grids is 70 × 70 × 60 (X × Y × Z) . The model domain initially consists of the upper crust with depths from 0 to 16 km, lower crust with depths from 16 to 32 km, and mantle with depths from 32 to 200 km. The upper crust is assumed to be a thermally conductive layer, while the lower crust and mantle are assumed to be convective regions. Different values for both the density and thermal conductivity are used for the thermally conductive layer and convective region. In the thermally conductive layer, the density and thermal conductivity are kept constant at 2600 kg/m 3 and 2.5 W/m•K, respectively. On the other hand, the values in the convective region are defined as functions of temperature and depth. Additionally, in the numerical simulation, a cold nose is not defined in the mantle wedge.
The initial temperature distribution is assumed to be represented by a half-space cooling model S1 with no flow. The boundary conditions for mantle flow are assumed to consist of no flow at the model surface (-z) and permeable at the five other boundaries (Fig. S5). The boundary condition for temperature is defined as a function of the depth and age of the oceanic plate at the trench in the vertical plane (+x), where the subduction of the oceanic plate begins S2 . The temperature of the model surface (-z) is held constant at 0°C, and the four other boundaries are set to be adiabatic. The geometry of the upper surface of the subducting plate was created by using the Slab2 model S3 . However, the results of the seismic velocity structure survey S4 showed that the depth of the plate surface was approximately 2 km shallower near the trench than that of the Slab2 model, so we set the slab geometry as such. Furthermore, in general, the iso-depth contours of the plate interface are aligned approximately parallel to the trench axis for most of the area but are perpendicular to the trench axis in the northeastern part of the model domain ( Fig. S1(a)). The locations of the hypocentres were used to create the slab geometry of Slab2. However, the number of earthquakes is small in the northeastern part of the model domain, so the reliability of the slab geometry may not be high ( Fig. S1(b)). Therefore, we changed the slab geometry in this region so that the isodepth contours were closer to the direction of the trench axis (Text S2 and Fig. S1(d)).
The thickness of the subducting plate was defined as a function of the age of the oceanic plate at the trench S5 . Using these criteria, we fixed the geometry of the top and bottom surfaces of the subducting oceanic plate and used this as a prescribed guide. We achieved the subduction of the oceanic plate by gradually pouring oceanic plate material into the guide at the convergence rates S6 . The convergence rates were determined by considering the subduction history based on the past plate rotation model S7 (Table S1).
The simulation period was set to 18 Myr, which begins from the start of plate subduction in the model domain to the present when the model domain, including the subducting plate, mantle wedge, and continental plate, reaches a nearly steady thermal state. To constrain the obtained thermal structure at the present (0 Ma), we used the observed heat flow data S8, S9 and constructed a model in which the residuals between the observations and calculations can be minimized ( Figure S6).

Effect of the subducted Yakutat terrane geometry on the results
The slab geometry used in Iwamoto et al. (2022) S7(b)). Therefore, the maximum dehydration gradient is approximately -2.0 wt%/km near the downdip in the southwestern part of the tectonic tremor-occurring area, whereas it is approximately -1.0 wt%/km near the centre of the tectonic tremor-generating area in the northeastern part ( Fig. S8(b)). We calculated the vertical sum of the dehydration gradient from the slab surface to the slab Moho and found that its spatial pattern is oblique to the tectonic tremor-occurring area and that the values are larger in the southwestern part of the tectonic tremor-generating area (Fig. S9). It is generally believed that the dehydrated water is brought to the plate interface and causes tectonic tremors, so it is difficult to explain the relationship between the two if their spatial distributions are oblique.
Therefore, we considered changing the slab geometry, which affects the thermal structure. In general, isodepth contours of the slab surface are approximately parallel to the trench axis, but as mentioned above, isodepth contours are oriented almost perpendicular to the trench axis in the northeastern part of the model (Fig. S1(a)). The hypocentre locations are used to determine the slab geometry of Slab2. In the northeastern part of the model, the number of earthquakes is small, so the reliability of the slab geometry may not be high ( Fig. S1(b)). First, we created a slab geometry (Fig. S1(c)) by modifying the slab geometry of Iwamoto et al. (2022) S10 ( Fig. S1(a)) so that the isodepth contours in the red dashed box can be almost parallel to the trench axis. The new slab geometry was then created by taking the weighted average of 8:2, 5:5, and 2:8 for the depth distributions of Fig. S1(a) and Fig. S1(c). Among these three slab geometries, we selected the most suitable one, for which the weighted RMS of the residuals between the observed and calculated heat flow was the smallest and the region with a large dehydration gradient was not oblique to the tectonic tremor occurring area. As a result, the geometry with a 2:8 weighted average satisfied these conditions ( Fig. S1(d)). In the Results and Discussion section of the main text, we showed the results of numerical simulations using this slab geometry.          Mapping Tools (GMT) S13 (version: GMT 4.5.7, URL link: https:// www. genericmapping-tools.org/download/).